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By means of variational methods and systematic numerical analysis, we demonstrate the existence 
of metastable solitons in three-dimensional (3D) free space, in the context of binary atomic conden¬ 
sates combining contact self-attraction and spin-orbit coupling, which can be engineered by available 
experimental techniques. Depending on the relative strength of the intra- and inter-component at¬ 
traction, the stable solitons feature a semi-vortex or mixed-mode structure. In spite of the fact that 
the local cubic self-attraction gives rise to the supercritical collapse in 3D, hence the setting produces 
no true ground state, the solitons are stable against small perturbations, motion, and collisions. 
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Introduction and model — Solitons result from the 
balance between dispersion and nonlinearity in diverse 
physical systems. Stable solitons in one dimension (ID) 
have been studied extensively in diverse media, most no¬ 
tably nonlinear optics and atomic Bose-Einstein conden¬ 
sates (BECs) [l|. Multidimensional solitons were also 
predicted to exist in ferromagnets 0, superconductors 
3, semiconductors Q, BECs baryonic matter Q, 
field theory Q, etc. However, creation of 2D and 31) 
bright solitons is a much more challenging problem than 
in ID. The fundamental difficulty is the fact that the 
ubiquitous cubic local self-attractive nonlinearity gives 
rise to the critical and supercritical collapse (blowup) 
in the 2D and 3D geometry, respectively |8l-ll0l|. which 
makes all the bright solitons unstable (the self-repulsive 
nonlinearity supports stable 2D dark solitons in the form 
of delocalized vortices EH)- Several theoretical schemes 
have been elaborated for the stabilization of 2D and 3D 
solitons. They rely on the use of trapping potentials 
(T^ - [l^ . sophisticated nonlinear interactions or 

nonlocal nonlinearity [^. However, it is commonly 
believed that a local cubic self-attraction may never give 
rise to stable solitons in 3D free space [1^ . 

Recently, an essential result [2^, which helps to re¬ 
solve a related but easier problem of the stabilization of 
solitons in 2D free space with local cubic attraction, has 
been reported in the framework of the model of a binary 
BEG subject to the action of spin-orbit coupling (SOC) 

a (solitons in ID SOC models have been predicted too 
, but their stability is obvious). It was found that the 
system gives rise to completely stable 2D bright solitons 
as the ground state (GS). The stabilization is explained 
by the fact that the linear SOC terms come with a co- 
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efficient whose dimension is inverse length. The usual 
2D systems without SOC feature a specific scaling invari¬ 
ance, which is closely related to the critical collapse. The 
scaling invariance makes the family of 2D solitons degen¬ 
erate (they are called Townes solitons in that case [27|). 
with a single value of the norm that does not depend on 
the soliton’s chemical potential. This norm determines 
the threshold for the onset of the critical collapse . 
Breaking the scaling invariance by introducing a fixed 
length scale leads to the stabilization of 2D solitons. This 
can be achieved by adding trapping potentials (l2l - [l^ or, 
in the free space, with the help of SOC [2^, which creates 
the missing GS by pushing the norm of the 2D solitons 
below the collapse threshold. A similar mechanism en¬ 
ables the stabilization of 2D spatiotemporal solitons in a 
planar optical coupler [2^, with the coupling’s temporal 
dispersion [2^ emulating the SOG effect. 


It has been previously shown that, besides the stabi¬ 
lization of 2D solitons, the interplay of SOC and intrinsic 
BEC nonlinearity give rise to a variety of other remark¬ 
able phenomena [30j. However, the possibility of sta¬ 
bilizing 3D solitons in free space with the help of SOC 
remained an open question. The fundamental difficulty 
is that, on the contrary to the 2D situation, the super¬ 
critical collapse in 3D has zero threshold, hence the norm 
cannot take values below the threshold, making the stabi¬ 
lization mechanism outlined above irrelevant in 3D. The 
present work reveals that, nevertheless, the self-attractive 
binary SOC condensate can support (meta)stable 3D 
solitons in free space, in spite of the fact that the set¬ 
ting has no GS at any value of the norm (in other words, 
the energy is unlimited from below). We find that the 
SOC-induced modffication of the dispersion of the 3D 
condensate may balance the attractive nonlinearity, cre¬ 
ating metastable solitons. In addition to the absence of 
the GS, another fundamental difference of this mecha¬ 
nism from what is outlined above for 2D is that the sta¬ 
bility of the 3D solitons is controlled not by the norm. 





2 


but rather by their energy. 

We follow the usual mean-field approach, defining 
dt(r) = {4’+, 4’-Y' the condensate wave function, with 
± referring to two pseudo-spin components. Fixing by 
means of rescaling the atomic mass and Planck’s con¬ 
stant to be 1 , we write the system’s energy as the sum of 
kinetic, SOC, and interaction terms: 


-^tot — -^kin d” ^soc d” -^int j (1) 

^kin = IJ £;soc=A j 

Eint = ~\ j {\^+\^ + \^-t + : 


where cr = {(Jx,<^y,<^z) are Pauli matrices, and p = — iV 
is the momentum operator. We adopt the 3D isotropic 
form of the SOC with strength A [^. The intra- and 
inter-component interaction strengths are defined, re¬ 
spectively, as —g and —rjg, with g > 0 corresponding to 
the self-attraction, 77 being the relative cross-nonlinearity 
strength. Below, we fix the nonlinearity strength, by 
rescaling the wave functions, to 5 = 1 and vary the SOC 
strength A, norm N, and cross-nonlinearity strength 77 . 

Dimensional analysis — If L is a characteristic size 
of the self-trapped condensate, an estimate for the 
amplitudes of the wave functions with norm N = 
J d^r (It/j+P- p |7/)_|2) is ~ There¬ 

fore, the three terms in Eq. © scale with L as 


Etot/N - CkinT-'-CsocAL-i-(c|^f^^ + 7VT-3 , 

( 2 ) 

with positive coefficients Ckin, Ooc, and ^ 

shown in Fig. [U Eq. ([2]) gives rise to a local minimum of 
Etot{L) at finite L, provided that 


0 < AfV < 



(cross) 

'"int 


v) 



( 3 ) 


Although this minimum cannot represent the GS (which 
formally corresponds to Etot —?► — 00 at L —>■ 0 in the 
collapsed state, i.e., the system has no true GS), it cor¬ 
responds to a self-trapped state stable against small per¬ 
turbations. Previously, a similar approximate analysis 
has correctly predicted stable quasi-2D solitons in dipo¬ 
lar BEC [ill. 

Condition suggests that metastable 3D solitons 
may exist in free space when the SOC term is present, 
while its strength A is not too large, N and 77 being not 
too large either. We confirm these expectations below by 
means of accurate numerical analysis. 

The Gross-Pitaevskii equation — Energy functional o 
gives rise to the Gross-Pitaevskii equation (CPE) for the 
spinor wave function. 
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FIG. 1: (Color online) Atot as a function of condensate’s size 
L, as per Eq. 0. The red solid, blue dashed, and green 
dot-dashed lines represent the energy’s variation when A = 0, 
and A > 0 does or does not satisfy condition ®, respectively. 


Assuming axial symmetry of the expected self-trapped 
states (it is the highest symmetry admitted by the SOC 
[^ 1 and using cylindrical coordinates {r,z,(j)), the sta¬ 
tionary wave function with integer vorticity m > 0 and 
chemical potential /i is looked for as 


'p- 




e™'^/i(r,z) A 

e*(™+i)'^/ 2 (r,z) ) ■ 


( 5 ) 


Following the terminology introduced for 2D solitons in 
Ref. [ill, self-trapped states (IS|) with m = Q are called 
semi-vortices (SVs), the states with ttt, > 1 being their 
excited states. Similar to the 2D system [2|, our calcula¬ 
tions demonstrate that the energy of the SV with m = 0 
is always lowest, therefore we focus on m = 0 . 

Due to the up-down symmetry of underlying Hamilto¬ 
nian o, degenerate to SV © is its flipped counterpart, 




e-*(™+i)'^/ 2 *(r,z) 
e"™'^ A*(r, z) 


( 6 ) 


with * standing for the complex conjugate. Although 
the system is axially symmetric, stationary states do not 
necessarily follow this symmetry. In particular, any su¬ 
perposition of ansdtze © and © breaks the symmetry. 
Following the nomenclature introduced in Ref. [2|, we 
call the state generated by such a superposition a mixed 
mode (MM). Approximating it by the superposition with 
mixing angle 9 [32j |. 


7 /)+ = (cos6») /i(r, z) - (sin6») f^{r, z) e , 

= (sin 6 ») f*{r, z) -b (cos 6 ») / 2 (r, z) e"^ , 


straightforward calculation relates its energy to that of 
the respective SV: 

Amm = Esv + (1 — ?7)sin^6lcos^0 AE, (8) 

AE = 2Trg J rdr (l/il'^-b 1 / 21 "^ - 4 |/ip|/ 2 p) . 

Our numerical calculations show that AE is always pos¬ 
itive, hence, like in the 2D case [2|, the SV (MM) has 
lower energy at 77 < 1 (77 > 1) . This prediction is con¬ 
firmed below by the full numerical analysis. 
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Variational analysis — To produce analytical results 
in a more accurate form than given by Eq. we here 
adopt the following ansatz for the SV: 

/„ = (A„ + iBr^z) (n = 1 , 2 ), 


branch, which does not satisfy the VK criterion, repre¬ 
sents solitons corresponding to the energy maximum on 
the blue dashed curve in Fig.[TJ In the limit of ^ — oo, 
they carry over into the well-known strongly unstable 3D 
solitons of the GPE [sj. 


with real parameters An, Bn, and > 0, /3„ > 0. The 
substitution of this ansatz into expression ([T|) for the full 
energy and minimizing it with respect to the free param¬ 
eters produces algebraic equations which can be readily 
solved numerically. Stable solitons correspond to finite 
values of an and /3„, while —>■ 0 (spreading) and 

an,l3n —>■ oo (collapsing) indicate that no solitons ex¬ 
ist. Results of the calculations are summarized in Fig. [21 
in which the stable 3D solitons are predicted to exist in 
the shaded areas. We thus conclude that the solitons in¬ 
deed exist, provided that A, N and rj are not too large, in 
agreement with the qualitative prediction of Eq. © from 
the dimensional analysis. In particular, an important 
conclusion is that, for fixed A and rj, the stable solitons 
always exist in a finite interval of the norm, 
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FIG. 3: (Color online) (a) Energies of the SVs and MMs, 
as predicted by the variational approach for p = A = 1 and 
N = 8. The two curves cross at t; = 1, where the SV and MM 
have equal energies, (b) The numerically (blue circles) and 
variationally (red squares) found chemical potential vs. the 
norm for the SVs at gr = A = 1 and y = 0.3. The numerical 
branch extends up to V = A^max, in agreement with Eq. ®. 


0<A^<A^max(A,?7) • (9) 

Furthermore, as shown in Fig.|3Ka), for 77 < I the energy 
of the SV predicted by the variational analysis (VA) is 
lower than that for the MM, and vice versa for 77 > I, in 
agreement with the prediction of Eq. ®. 



FIG. 2: (Color online) 3D stable solitons are predicted by the 
variational calculation in blue shaded regions of the respective 
parameter planes. In panel (a), these are SVs (semi-vortices) 
at 77 < 1, and MMs (mixed modes) at 77 > 1, with the bound¬ 
ary between them depicted by the black solid line. In (b), the 
entire stability area is filled by the solitons of both types, as 
the SVs and MMs have equal energies at 77 = 1. The predic¬ 
tions are accurately confirmed by full numerical simulations, 
as indicated by red crosses and black dots, which indicate, 
respectively, the absence and presence of stable solitons for 
respective sets of parameters. 


The red squares in Fig. |2Kb) represent the variational 
results for the soliton’s chemical potential, /i, plotted 
as a function of norm N for g = X = 1 and 77 = 0.3. 
In agreement with the analytical prediction given by 
Eq. (|21), there is no threshold (minimum norm) neces¬ 
sary for the appearance of the solitons, which exist up to 
a, N = A^max- Furthermore, the negative slope of the de¬ 
pendence, dg/dN < 0, of the upper branch is an indica¬ 
tion of the stability of the soliton families, pur suant to the 
Vakhitov-Kolokolov (VK) criterion [sl. [^.13^ . The lower 


Full numerical calculations — The prediction for the 
existence of the stable 3D solitons in free space, provided 
by the analytical approximations, calls for verification 
by direct simulations of GPE ©. First, we generated 
stationary states by running the simulations in imagi¬ 
nary time. Typical examples of the so produced SV and 
MM density profiles are displayed in Fig. U) Symbols in 
Fig.m which indicate the absence and presence of stable 
solitons, are in good agreement with the VA. 



-3 -3 -3 -3 


FIG. 4: (Golor online) Density profiles of 3D solitons for N = 
8 and g = X = 1. (a) An SV for 77 = 0.3, whose fundamental 
and vortical components, \ip+\ and \ip-\, are plotted in (al) 
and (a2), respectively, (b) An MM for 77 = 1.5, with (bl), 
(b2) displaying \ip+\ and \ip-\, respectively. In each subplot, 
different colors represent constant-magnitude surfaces, \ip± \ = 
(0.96,0.4,0.04) X |7/7±|,n,x. 

The blue circles in Fig. EKb) represent the numerically 
obtained chemical potential, which are in good agreement 
with the prediction of the VA. The unstable branch from 
the VA, however, cannot be produced by the imaginary- 
time integration. We have verified the stability of the 
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solitons belonging to the upper branch in Fig. IH^b) by 
real-time simulations with random perturbations added 
to the initial conditions, confirming that the VA accu¬ 
rately predicts the SV and MM stability areas which are 
displayed in Figs. [ 2 ] and [S] 

Setting quiescent solitons in motion is another nontriv¬ 
ial issue, as the SOC terms break the Galilean invariance 
of the system. To construct solitons moving along the z- 
axis with velocity Uj, so that '0± = 4>, z — Vzt, t), we 

have rewritten the GPE system (|T]) in the respective mov¬ 
ing reference frame. In this form, the velocity term affects 
the SOG strength along the z axis, breaking the symme¬ 
try between the two components of the spinor. As a re¬ 
sult, positive (negative) Vz tends to increase the popula¬ 
tion of the spin-down (-up) component. In Fig.[5l we plot 
the ratio of the spin populations as a function of Vz ■ Both 
VA and numerical results are displayed, showing qualita¬ 
tively similar results. At Vz < —0.9 and Vz > -1-0.4, the 
moving semi-vortex practically degenerates into a single¬ 
component soliton - the fundamental or vortical one, re¬ 
spectively - thus reducing the setting to that for the sin¬ 
gle GPE with the cubic self-attraction, where all 3D soli¬ 
tons are strongly unstable. Consequently, the speed of 
the stably moving solitons cannot be too large. 



FIG. 5: (Color online) The ratio of the spin populations as 
a function of velocity Vz for the moving SV with N = 8, 
g — X = 1, rj — 0.3, and N± = J |i/)±(r)p. The red 
dashed lines with squares are variational results, while the 
blue solid lines with circles are obtained numerically, using 
the imaginary-time integration in the moving reference frame. 

Finally, to consider collisions between moving solitons, 
we place two solitons centered at initial positions (r, z) = 
( 0 , ±Zo)j and include a trapping potential, f 2 ^(r^ -|-z^)/ 2 . 
The solitons then start moving to collide at the trap cen¬ 
ter, with the trapping frequency D used to control the 
collision velocity. Figure [S] depicts two collision events 
for the same initial soliton pair. In panel (a), the slowly 
moving solitons feature a quasi-elastic collision, while in 
(b) the collision leads to destruction of faster solitons. 
This shows the solitons are robust against slow collisions. 

Conclusion — The combination of the analytical and 
numerical methods reveals that stable free-space 3D soli¬ 
tons can be supported in the binary atomic conden¬ 
sate with attractive interactions and properly engineered 
SOC, notwithstanding the presence of the supercritical 
collapse in the same setting. This is the first example of 
metastable solitons in the 3D homogeneous environment 
with local cubic self-attraction, which exist in spite of the 
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FIG. 6: (Golor online) Gollisions of stable 3D SVs in the 
harmonic trap for N = 8, g = X = 1, and rj = 0.3. Panels 
(al, a2), (a3, a4), (a5, a6), (a7, a8) and (a9, alO) display 
density distributions for Q = 0.5 at t = 0, 1.2, 3.2, 4 and 6, 
respectively. Panels (bl, b2), (b3, b4), (b5, b6), (b7, b8) and 
(b9, blO) display the distributions for D = 1 at t = 0, 1, 1.4, 
2 and 3, respectively. In all panels, the left and right subplots 
display, severally, \ip+\ S'Hd iV'-l- 


nonexistence of the GS (ground state) in the system. The 
SOC plays a crucial role for the stabilization, altering the 
energy of the self-trapped states so as to create the lo¬ 
cal energy minimum. This is the fundamental difference 
from the recently discovered stabilization mechanism in 
2D [^, which readily creates a missing GS below the 
critical value of the norm (at N < N^), where solitons, 
if any, cannot be destabilized by the critical collapse, as it 
does not occur at V < Na, but no solitons could be cre¬ 
ated at V > Ncr- In 3D; the existence of the metastable 
solitons is controlled not by the norm [in an appropriate 
parameter region, they can be created for any TV, al¬ 
though the appropriate region becomes very narrow for 
very large TV, as seen in Fig. |2Ib)], but by the energy, as 
the above analysis clearly shows. 

Although we have adopted the isotropic SOC term in 
the Hamiltonian, in the form of Ap • cr, the stabilization 
of the 3D solitons does not critically depend on this form, 
additional analysis demonstrating that the metastable 
3D solitons exist as well if the SOC strength is different 
along different axes. It may also be interesting to find out 
if 3D solitons can be stabilized by spatially localized SOC 
(for ID solitons, this setting was studied in Ref. [s^, but 
the stability is not an issue in that case). Influence of the 
Zeeman splitting, which breaks the up-down symmetry 
of the spinor components, on the stability of the solitons 
is another relevant problem for further analysis. 

On the experimental side, 2D SOC was recently cre¬ 
ated in an ultracold Fermi gas [s^. Realization of 3D 
SOC may be expected in the near future, as there is no 
fundamental obstacle for doing that. 
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